{ "cells": [ { "cell_type": "markdown", "id": "cf5a3b29", "metadata": {}, "source": [ "# User Defined Reactions in MUSICA\n", "There is an additional type of reaction in MUSICA called User Defined reactions, which are a more flexible type of reaction that supports defining per-grid-cell reaction rates.\n", "This differs from the traditional reaction types where the reaction rates are calculated.\n", "Defining the reaction rate can be useful in cases where you have observations and need to find the rate that most closely matches those observations.\n", "There are multiple classes of User Defined reactions, with there being the default UserDefined class as well as the Emission and FirstOrderLoss classes also being User Defined under the hood." ] }, { "cell_type": "markdown", "id": "63424343", "metadata": {}, "source": [ "## 1. Importing Libraries\n", "Below is a list of the required libraries for this tutorial:" ] }, { "cell_type": "code", "execution_count": 1, "id": "f81fdb7f", "metadata": {}, "outputs": [], "source": [ "import musica\n", "import musica.mechanism_configuration as mc\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import qmc\n", "import pandas as pd\n", "import numpy as np\n", "import seaborn as sns\n", "pd.set_option('display.float_format', str) # This is done to make the arrays more readable\n", "np.set_printoptions(suppress=True) # This is done to make the arrays more readable" ] }, { "cell_type": "markdown", "id": "98b68587", "metadata": {}, "source": [ "## 2. Redefining the Original Species and Reactions\n", "As with the previous two tutorials, this is simply defining the three species and two reactions that have been used already." ] }, { "cell_type": "code", "execution_count": 2, "id": "50d98bb3", "metadata": {}, "outputs": [], "source": [ "A = mc.Species(name=\"A\")\n", "B = mc.Species(name=\"B\")\n", "C = mc.Species(name=\"C\")\n", "species = [A, B, C]\n", "gas = mc.Phase(name=\"gas\", species=species)\n", "\n", "r1 = mc.Arrhenius(\n", " name=\"A_to_B\",\n", " A=4.0e-3, # Pre-exponential factor\n", " C=50, # Activation energy (units assumed to be K)\n", " reactants=[A],\n", " products=[B],\n", " gas_phase=gas\n", ")\n", "\n", "r2 = mc.Arrhenius(\n", " name=\"B_to_C\",\n", " A=4.0e-3,\n", " C=50, \n", " reactants=[B],\n", " products=[C],\n", " gas_phase=gas\n", ")" ] }, { "cell_type": "markdown", "id": "d06d4447", "metadata": {}, "source": [ "## 3. Defining User Defined Reactions\n", "This code cell defines three new reactions: one User Defined, one Emission, and one First Order Loss.\n", "It then creates a new mechanism containing both the old reactions as well as the three new reactions." ] }, { "cell_type": "code", "execution_count": 3, "id": "6781a5a5", "metadata": {}, "outputs": [], "source": [ "r3 = mc.UserDefined(\n", " name=\"complex_rxn\",\n", " scaling_factor=1.0,\n", " reactants=[A, (2, B)],\n", " products=[A, (2, C), B],\n", " gas_phase=gas\n", ")\n", "\n", "r4 = mc.Emission(\n", " name=\"B_emission\",\n", " scaling_factor=1.0,\n", " products=[B],\n", " gas_phase=gas\n", ")\n", "\n", "r5 = mc.FirstOrderLoss(\n", " name=\"C_loss\",\n", " scaling_factor=1.0,\n", " reactants=[C],\n", " gas_phase=gas\n", ")\n", "\n", "mechanism = mc.Mechanism(\n", " name=\"musica_micm_example\",\n", " species=species,\n", " phases=[gas],\n", " reactions=[r1, r2, r3, r4, r5]\n", ")" ] }, { "cell_type": "markdown", "id": "da99a490", "metadata": {}, "source": [ "## 4. Creating the Solver, State, and Latin Hypercube Sampler\n", "This code is a rehash that includes code similar to the previous Hypercube tutorial, just with the LHS being expanded to 8 dimensions to account for the three new User Defined reactions as the last three dimensions. " ] }, { "cell_type": "code", "execution_count": 4, "id": "3947b2c3", "metadata": {}, "outputs": [], "source": [ "solver = musica.MICM(mechanism = mechanism, solver_type = musica.SolverType.rosenbrock_standard_order)\n", "\n", "num_grid_cells = 100\n", "state = solver.create_state(num_grid_cells)\n", "\n", "ndim = 8\n", "nsamples = num_grid_cells\n", "\n", "# Create a Latin Hypercube sampler in the unit hypercube\n", "sampler = qmc.LatinHypercube(d=ndim)\n", "\n", "# Generate samples\n", "sample = sampler.random(n=nsamples)\n", "\n", "# Define bounds for each dimension\n", "l_bounds = [275, 100753.3, 0, 0, 0, 0, 0, 0] # Lower bounds\n", "u_bounds = [325, 101753.3, 20, 10, 20, 1, 0.5, 0.5] # Upper bounds\n", "\n", "# Scale the samples to the defined bounds\n", "sample_scaled = qmc.scale(sample, l_bounds, u_bounds)\n", "\n", "temperatures = sample_scaled[:, 0]\n", "pressures = sample_scaled[:, 1]\n", "concentrations = {\n", " \"A\": [],\n", " \"B\": [],\n", " \"C\": []\n", "}\n", "concentrations[\"A\"] = sample_scaled[:, 2]\n", "concentrations[\"B\"] = sample_scaled[:, 3]\n", "concentrations[\"C\"] = sample_scaled[:, 4]" ] }, { "cell_type": "markdown", "id": "3af454f9", "metadata": {}, "source": [ "## 5. Creating and Setting the User Defined Parameters\n", "Next, each of the new reactions is given a dictionary with a name as the key and their respective LHS outputs as the value.\n", "For these reactions, the name must match the name given to the reactions above, with an added prefix, being:\n", "* \"USER.\" for User Defined reactions,\n", "* \"EMIS.\" for Emission reactions, and\n", "* \"LOSS.\" for First Order Loss reactions.\n", "\n", "After that, the state's set_user_defined_rate_parameters() function is called to add all of the reaction rates to the system." ] }, { "cell_type": "code", "execution_count": 5, "id": "49b15078", "metadata": {}, "outputs": [], "source": [ "complex_rates = {\"USER.complex_rxn\": sample_scaled[:, 5]}\n", "emission_rates = {\"EMIS.B_emission\": sample_scaled[:, 6]}\n", "loss_rates = {\"LOSS.C_loss\": sample_scaled[:, 7]}\n", "\n", "state.set_user_defined_rate_parameters(complex_rates)\n", "state.set_user_defined_rate_parameters(emission_rates)\n", "state.set_user_defined_rate_parameters(loss_rates)" ] }, { "cell_type": "markdown", "id": "d746c84f", "metadata": {}, "source": [ "## 6. Setting the Conditions, Running the Solver, and Visualizing the Data\n", "As with step 4, this is a copy of some previous steps from the Hypercube tutorial, with the only change being the simulation length being extended to 100 seconds.\n", "Note the extreme difference in the concentration curves after the three new reactions were added; B is now the dominant species instead of C after approximately 50 seconds of the simulation taking place." ] }, { "cell_type": "code", "execution_count": 6, "id": "dda25f0a", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | time.s | \n", "ENV.temperature.K | \n", "ENV.pressure.Pa | \n", "ENV.air number density.mol m-3 | \n", "CONC.A.mol m-3 | \n", "CONC.B.mol m-3 | \n", "CONC.C.mol m-3 | \n", "
---|---|---|---|---|---|---|---|
0 | \n", "0 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "2.893750993982109 | \n", "2.673213216234265 | \n", "9.433830407714845 | \n", "
1 | \n", "0 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "7.053727382338282 | \n", "2.260055721965611 | \n", "4.4571713767278744 | \n", "
2 | \n", "0 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "10.711091077403818 | \n", "6.697825729904961 | \n", "10.843190900041073 | \n", "
3 | \n", "0 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "15.546051530784109 | \n", "7.64337289065524 | \n", "1.209212898549111 | \n", "
4 | \n", "0 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "14.515576668606005 | \n", "6.473953665388967 | \n", "9.338659341590636 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
10095 | \n", "100 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "5.934136681910782e-10 | \n", "78.38313350457263 | \n", "0.9507647383147676 | \n", "
10096 | \n", "100 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "3.5536302092439606e-10 | \n", "93.3521115434026 | \n", "0.9755165898633137 | \n", "
10097 | \n", "100 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "6.772157344069222e-10 | \n", "23.47690407300151 | \n", "0.6340441771972303 | \n", "
10098 | \n", "100 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "4.7929670861011336e-11 | \n", "37.75925021640928 | \n", "0.849375775534841 | \n", "
10099 | \n", "100 | \n", "309.7352766001884 | \n", "101389.99779492112 | \n", "39.37043899688805 | \n", "3.5291165084174213e-10 | \n", "101.85877228753668 | \n", "3.3624463280973975 | \n", "
10100 rows × 7 columns
\n", "